function [A_hat,Lambda_hat,logLik,BIC,Q] = MLE_block_cov(X,vB)

% Inputs:
% X - T x n data array of n-dimensional vectors, where T is a sample size .
% It is assumed that Xt is i.i.d. N(0, Sigma), t = 1,...,T, and Sigma has K x K block structure.
% vB - K-dimensional vector of block sizes, where K is a number of blocks (clusters).
% Entries of vB must be positive integers and sum up to n.
% Note: data in matrix X must be pre-ordered in such a way that variables
% in the first vB(1) columns correspond to the first cluster,
% variables in the next vB(2) columns are from the second cluster, etc .

% Output:
% A_hat - ML estimate of K x K matrix A from the canonical representation
% Lambda_hat - K-dimensional vector of ML estimates of lambda_1,...,lambda_K
% logLik - value of log-likelihood function at optimum
% BIC - Bayesian information criterion
% Q - n x n orthogonal matrix from the canonical representation

[T,n] = size(X); K = length(vB);

% Check if input is of proper format : vB consists of positive integers and sum(vB) equals n
if all(vB > 0) && all(mod(vB,1) == 0) && (sum(vB) == n)
    
    vB_ind = [0; cumsum(vB)];   % Get indicies of the block partition implied by vB
    Q = zeros(n);               % Allocate n � n matrix for Q

    % Estimate A
    Y0 = zeros(T,K);
    for k = 1:K
        % Construct k-th column of Q
        Qk = zeros(n,1); vk = ones(vB(k),1)/sqrt(vB(k));
        Qk(vB_ind(k)+1:vB_ind(k+1)) = vk;

        Q(:,k) = Qk;      % Update Q
        Y0(:,k) = X*Qk;   % Transform data vectors
    end
    A_hat = Y0'*Y0/T; % MLE of A

    % Estimate lambda_1,...,lambda_K
    Lambda_hat = zeros(K,1);
    for k = 1:K
        % Explicitly get an orthogonal complement to vk and restore remaining columns of Q
        Qk_ort = zeros(n,vB(k)-1); vBk_lin = (0:vB(k))';
        vk_ort = triu(ones(vB(k))*diag([0;1./sqrt(vBk_lin(2:end-1).*vBk_lin(3:end))]));
        vk_ort(eye(vB(k),'logical')) = -diag(vk_ort).*vBk_lin(1:end-1);
        Qk_ort(vB_ind(k)+1:vB_ind(k+1),:) = vk_ort(:,2:end);

        Q(:,K-k+vB_ind(k)+2:K-k+vB_ind(k+1)) = Qk_ort;   % Update Q
        Yk = X*Qk_ort;                                   % Transform data vectors
        Lambda_hat(k) = trace(Yk*Yk')/((vB(k)-1)*T);     % MLE of lambda_k
    end
    
    % Compute log-likelihood and BIC
    logLik = -0.5*n*T*(log(2*pi)+1)-0.5*T*(log(det(A_hat))+nansum((vB-1).*log(Lambda_hat)));
    BIC    = (0.5*K*(K+1)-sum(vB==1))*log(n*T)-2*logLik;

else
    fprintf ('Error : input is of wrong format \n'); return
end